/*
 * ReadBoundaryConditionforThermalResponse.C
 *
 *  Created on: Jan 11, 2020
 *      Author: liuxiao
 */

#include "libmesh/system.h"
#include "libmesh/radial_basis_interpolation.h"
#include<iostream>

#include "../../include/materials/ReadBoundaryConditionforThermalResponse.h"
#include "libmesh/quadrature.h"
using std::cout;
using std::endl;
registerMooseObject("TrilobitaApp", ReadBoundaryConditionforThermalResponse);
template<>
InputParameters validParams<ReadBoundaryConditionforThermalResponse>()
{
  InputParameters params = validParams<Material>();

  params.addRequiredCoupledVar("wall_temperature_variable", "Coupled Wall Temperature");
  params.addParam<Real>("Pe_O2", 0.21, "partial pressure of O2 at boundary layer edge,1");
  params.addParam<Real>("Pe_N2", 0.79, "partial pressure of N2 at boundary layer edge,1");
  params.addParam<Real>("M_e",0.029, " average molar mass of gas at boundary layer edge,kg/mol");
  params.addParam<std::string>("flux_data", "flux data from file");
  params.addParam<Real>("power", 2, "distance power");
  params.addParam<int>("interp_pts", 8, "number of interpolate points");
  params.addParam<std::vector<Real> >("scale",  "scale");
  params.addParam<Real>("fluxcoff",  "fluxcoff");
  return params;
}

ReadBoundaryConditionforThermalResponse::ReadBoundaryConditionforThermalResponse(const InputParameters & parameters) :
      Material(parameters),
	  _qc(declareProperty<Real>("coldwall_heatflux")),
	  _qc_0K(declareProperty<Real>("heatflux_0K")),
	  _hr(declareProperty<Real>("recovery_enthalpy")),
	  _Pe(declareProperty<Real>("boundary_layer_edge_pressure")),
	  _Pe_atm(declareProperty<Real>("Pe_atm")),
	  _Pe_O2(declareProperty<Real>("Pe_O2")),
	  _Pe_N2(declareProperty<Real>("Pe_N2")),
	  _M_e(declareProperty<Real>("M_e")),
	  _Tw_qc(declareProperty<Real>("Tw_for_coldwall_heatflux")),
	  _Tw(coupledValue("wall_temperature_variable")),
	  _Pe_O2_p(getParam<Real>("Pe_O2")),
	  _Pe_N2_p(getParam<Real>("Pe_N2")),
	  _M_e_p(getParam<Real>("M_e")),
	  _flux_file(isParamValid("flux_data") ? getParam<std::string>("flux_data") : ""),
	  _scale(getParam<std::vector<Real> >("scale")),
	  _fluxcoff(getParam<Real>("fluxcoff"))
{
	readFile();

	Real power = getParam<Real>("power");
	int n_interp_pts = getParam<int>("interp_pts");
	_field_name.push_back("qc");
	_field_name.push_back("hr");
	_field_name.push_back("tw0");
	_field_name.push_back("pe");

	for(int t = 0  ; t < _num_time_step; ++t)
	{
		std::vector<Real> src_vals;
		for(int i = 0; i < _num_pts; ++i)
		{
			src_vals.push_back(_src_qc[t][i]*_fluxcoff);
			src_vals.push_back(_src_hr[t][i]);
			src_vals.push_back(_src_tw0[t][i]);
			src_vals.push_back(_src_pe[t][i]);
		}

		_idis.push_back(new InverseDistanceInterpolation<LIBMESH_DIM>(this->comm(), n_interp_pts, power));
		_idis[t]->set_field_variables(_field_name);
		_idis[t]->get_source_points() = _src_pts[t];
		_idis[t]->get_source_vals() = src_vals;
	}
}

void ReadBoundaryConditionforThermalResponse::readFile()
{
	using namespace std;
	ifstream qc_read(_flux_file.c_str());
    if(!qc_read.good())
	     mooseError("Error opening file '" + _flux_file + "' from  data.");

    string line;
    getline(qc_read, line);
    istringstream iss(line);
    Real f;
    vector<Real> head;
    while(iss >> f)
       head.push_back(f);

    _num_pts = head[0];
    _num_time_step = head[1];
    _src_qc.resize(_num_time_step);
    _src_hr.resize(_num_time_step);
    _src_tw0.resize(_num_time_step);
    _src_pe.resize(_num_time_step);
	_src_pts.resize(_num_time_step);
    for (int t = 0; t < _num_time_step; ++t)
    {
    	vector<Real> store;
    	for(int i = 0; i < _num_pts ; ++i)
    	{
    		store.resize(_num_pts);
    	   	int j = 0;
    		getline(qc_read, line);
    		istringstream transline(line);
    		Real d;
    		while (transline >> d)
    		{
    			store.push_back(d);
    			store[j] = d;
    			j+=1;
    		}

    	    Point pp(store[0]*_scale[0],store[1]*_scale[1],store[2]*_scale[2]);
    	    _src_hr[t].push_back(store[3]);
    	    _src_qc[t].push_back(store[4]);
    	    _src_tw0[t].push_back(store[5]);
    	    _src_pe[t].push_back(store[6]);
    	    _src_pts[t].push_back(pp);
    	}
    	_time_step.push_back(store[7]);
	}

}


void ReadBoundaryConditionforThermalResponse::interpolate(std::vector<Real> &qc, std::vector<Real>  &hr,std::vector<Real>  &tw0, std::vector<Real>  &pe,std::vector<Point> &pts, Real t)
{
	int lower, upper;
	if(t < _time_step.front())
	{
		lower = 0;
		upper = 0;
	}
	else if(t > _time_step.back())
	{
		lower = _num_time_step - 1;
		upper = _num_time_step - 1;
	}
	else
	{
	    std::vector<double>::iterator up = std::lower_bound(_time_step.begin(), _time_step.end(), t);
	    upper = static_cast<unsigned int>(std::distance(_time_step.begin(), up));
	    if (_time_step[upper] == t)
	      lower = upper;
	    else
	      lower = upper - 1;
	}
	Real lam = 0.5;
	if(lower != upper)
	{
		lam = (_time_step[upper] - t)/(_time_step[upper] - _time_step[lower]);
	}
	std::vector<Real> tgt_val_lower,  tgt_val_upper;
	_idis[lower]->interpolate_field_data(_field_name, pts, tgt_val_lower);
	_idis[upper]->interpolate_field_data(_field_name, pts, tgt_val_upper);

	std::vector<Real> tgt_val_current;
	for(unsigned int i  = 0; i < tgt_val_lower.size(); ++i)
	{
		tgt_val_current.push_back(lam*tgt_val_lower[i] + (1-lam)*tgt_val_upper[i]);
	}

	int num_val = _field_name.size();
	for(unsigned int ip = 0; ip < pts.size(); ++ip)
	{
		qc.push_back(tgt_val_current[num_val*ip]);
		hr.push_back(tgt_val_current[num_val*ip+1]);
		tw0.push_back(tgt_val_current[num_val*ip+2]);
		pe.push_back(tgt_val_current[num_val*ip+3]);

	}

}

void ReadBoundaryConditionforThermalResponse::computeProperties()
{
	std::vector<Real> qc,hr,tw_qc,pe;
	std::vector<Point> pts;
	for(unsigned int ip = 0; ip < _q_point.size(); ++ip)
	{
		pts.push_back(_q_point[ip]);
	}
	interpolate(qc,hr,tw_qc,pe,pts,_t);

	for(_qp = 0; _qp < _qrule->n_points(); ++_qp)
	{
		_qc[_qp]=qc[_qp];
		_Tw_qc[_qp]=tw_qc[_qp];
		_hr[_qp]=hr[_qp]*1000;
//		_qc_0K[_qp]=_qc[_qp]*_hr[_qp]/(_hr[_qp]-1000*_Tw_qc[_qp]);
		_qc_0K[_qp]=_qc[_qp];
		_Pe[_qp]=pe[_qp];
		_Pe_atm[_qp]=pe[_qp]/101325;
        _Pe_O2[_qp]=_Pe_atm[_qp]*_Pe_O2_p;
        _Pe_N2[_qp]=_Pe_atm[_qp]*_Pe_N2_p;
        _M_e[_qp]=_M_e_p;
	}
}



